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ABSTRACT 

The Rossiter-McLaughlin (hereafter RM) effect is a key tool for measuring the projected spin-orbit angle between stellar spin axes 
and orbits of transiting planets. However, the measured radial velocity (RV) anomalies produced by this effect are not intrinsic and 
depend on both instrumental resolution and data reduction routines. Using inappropriate formulas to model the RM effect introduces 
biases, at least in the projected velocity V sin compared to the spectroscopic value. Currently, only the iodine cell technique has been 
modeled, which corresponds to observations done by, e.g., the HIRES spectrograph of the Keck telescope. In this paper, we provide a 
simple expression of the RM effect specially designed to model observations done by the Gaussian fit of a cross-correlation function 
(CCF) as in the routines performed by the HARPS team. We derived also a new analytical formulation of the RV anomaly associated 
to the iodine cell technique. For both formulas, we modeled the subplanet mean velocity v p and dispersion /3 P accurately taking the 
rotational broadening on the subplanet profile into account. We compare our formulas adapted to the CCF technique with simulated 
data generated with the numerical software SOAP-T and find good agreement up to Vsinj* < 20 km.s -1 . In contrast, the analytical 
models simulating the two different observation techniques can disagree by about 10<x in Vsin/* for large spin-orbit misalignments. 
It is thus important to apply the adapted model when fitting data. 

Key words. Astronomical instrumentation, methods and techniques - Instrumentation: spectrographs - Methods: analytical - 
Methods: data analysis - Methods: numerical - Techniques: spectroscopic - Planetary systems - Planets and satellites: fundamental 
parameters - Planets and satellites: general - Stars: planetary systems 



1. Introduction 

Transiting planets produce radial velocity (RV) anomalies when 
crossing the disk of their star. This mechanism, known as the 
Rossiter-McLaughlin effect (hereafter RM effect), is due to the 
stellar proper rotation and the fact that during a transit, a planet 
successively covers different portions of the stellar disk with 
different average velocities along the line of sight ( |Holt| [1893 ; 

~ [T9l4] l. The RM effect has gained 



Rossiter 1924 McLaughlin 



importance in the exoplanet community since it allows the mea- 
surement of the projected angle between the stellar spin-axis and 
the orbit of the planet. 

The first measurements of the RM effect induced by a transit- 
ing planet were performed almost simultaneously with two dif- 
ferent instruments on the same bright star HD209458. |Queloz| 
et al. (2000) observed the signal with the ELODIE spectrograph 
on the 193cm telescope of the Observatoire de Haute Provence, 
while [Bundy & Marcy| p000) used the HIRES spectrograph on 
the Keck telescope. As explained later on, the choice of the in- 
strument and, more particularly, the subsequent reduction anal- 
ysis, have a non negligible impact on the resulting shape of sig- 
nal. It is thus interesting to see that the two kinds of instrument, 
coupled with their own data treatment, which are still employed 
today, have been used in parallel since the beginning. 



The first planet-host stars studied with this technique were 
all compatible with low obliquity. It was along the lines of the 
model of planet migration in a protoplanetary disk. But then, 
several misaligned systems have been detected, starting with 
XO-3 with an angle initially announced at 70° ± 15° (Hebrard 
[etaLl [2008) 1 and then refined to 37.3° + 3° fvVinn et al.| |2009l >7 
yet still significantly misaligned. With the growth of the sample, 
a first correlation appeared showing that the hottest stars tend 
to be more misaligned (Wi nn et"aL] |2010] l. Additionally, a first 
statistical comparison between observations and theoretical pre- 



dictions have been performed (Triaud et al. 2010 1, suggesting 



that some of the hot Jupiters might be the result of the interac- 
tion with a stellar companion, leading to a Lidov-Kozai mech- 
anism, characterized by phases of large eccentricity and incli- 
nation and followed by a circularization by tides raised on the 
planet as it approaches the star (|Wu & Murray| [2003 ; Fabrycky 



& Tremaine 2007| l. Other scenarios have then been developed 



to explain the formation of hot Jupiters, such as plane t -plane t 
scattering (Rasio & Ford 1996 Beauge & Nesvorny| 2012) , 
the crossing of secular resonances (Wu & Lithwick 2011 1, or 



the Lidov-Kozai mechanism produced by a planetary compan- 
ion ( |Naoz et aT] |2011[ |Nagasawa & Ida| |2011| ). A new trend 
between age and obliquity has also been found, suggesting that 
tidal dissipation may play an important role in the evolution of 
those systems (Triaud 2011) . 
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Accurate modelings of the RM effect are thus needed to get 
reliable information on the current obliquity of stars and to test 
theoretical predictions. In the literature, one can find several ana- 
lytical expressions to model this effe ct fKopal||1942||Ohta et al.| 

~ 2010| |2011) . They are not 



2005 ; Gimenez 2006 Hirano et al. 



all identical because they model different techniques of radial 
velocity measurements. 

This raises an issue that should be considered with caution. 
RVs measured by different techniques or even by different in- 
struments using the same algorithm can differ by more than a 
constant offset. To illustrate this point, we consider a random 
variable with a given probability distribution function (PDF), 
and ask what its average value is. The term average is vague 
and can mean different quantities: mean, median, mode, etc. 
Nevertheless, if the PDF is symmetrical, any of these quantities 
lead to the same result, eventually with different robustnesses 
with respect to noise. But if the PDF is not symmetrical, each 
estimator of the average value provides different results that can- 
not be compared directly. The situation is similar in RV obser- 
vations, as noticed by Hirano et al. (2010). There are at least 
two different ways to measure RVs. One relies on the iodine cell 
technique which consists in fitting an observed spectrum with 
a modeled one that is Doppler-shifted (Butler et al. 1996 ). The 



other is based on a Gaussian fit to a cross-correlation function 
(CCF) ( |Baranne et aT| [1996} [Pipe et al.||2002) . The former tech- 
nique is applied to observations made with HIRES on the Keck 
telescope or with HDS at the Subaru telescope, while the lat- 
ter is the routine of, e.g., SOPHIE at the Observatoire de Haute 
Provence or HARPS at La Silla Observatory. If stars are affected 
by neither spots nor transiting planets, their spectral lines have 
constant shapes, and thus the RVs derived by any instruments 
may only differ by a constant offset. In contrast, spectral lines 
are deformed during transits, and these deformations vary with 
time. As a consequence, each analysis routine is expected to lead 
to a different signal. In turn, it is important to have an analytical 
model adapted to each analysis routine in order to interpret RM 
data. Moreover, it also means that one should not combine RM 
measurements from different instruments. 

The first analytical expressions of the RM effect were de- 
rived by |Kopal| ( f!9"42"| >, |Ohta et al.| ( |2005) , and |Gimene^ ( |2l)06] l. 
They all computed the weighted mean velocity, hereafter called 
fmean, along the line of sight of the stellar surface uncovered by 
the planet. This mean velocity, weighted by the surface intensity, 
leads to an exact expression of the form 



/ 



1-/ 



(1) 



where / is the fraction of the flux blocked by the planet disk and 
Up is the average velocity of the surface of the star covered by the 
planet. This expression is simple and exact. There is no assump- 
tion behind it. But, it does not correspond to the quantity that 
is actually measured in the observations by either the iodine cell 
technique, or by the Gaussian fit of the CCF. Thus, the analytical 
prediction i> mea n is systematically biased when compared directly 
with observations. The difference increases for high stellar rota- 
tional velocities. 

To solve this problem, Hirano et al.| ( 2010| > propose a new 
analytical expression closely related to the reduction algorithm 
of the iodine cell technique (see Fig. [TJ). This method consists in 
fitting a shifted spectrum outside of transits, which is modeled, 
here, by a single averaged spectral line and notec^j (v - v), 

1 Here, and throughout the paper, the line profiles are expressed as a 
function of radial velocity v, instead of wavelength. 



with the line profile detected during a transit !Ftransit(fl)> via the 
Doppler shift v. The quantity provided by the formulas of Hirano 
|et al.| ( |2010| l is the value, hereafter denoted vhw, of the veloc- 
ity v that maximizes the cross-correlation between the two line 



profiles C(v) = J TstaAv - v) ( F\ IW sii(v) dv. We show in Sect. 3.1 
that if Tx m is an even function, the maximization of the cross- 
correlation C(v) is indeed identical to the minimization of the 
chi-square associated to the fit of Ttrmsit(v) by T s \m{v - !>)■ In that 
case, the result depends on the actual shape of the line profiles, 
as observed in practice. In the simplest case where the line pro- 
files of both the nonrotating and the rotating star are Gaussian, 
this method leads to ( |Hirano et al.||2010] l 



3/2 



fVp 



m +/%)) 



(2) 



where /?+ and /3 P are the dispersion of the Gaussian profiles 
Fstar(v) and Tpi a (v-v p ) respectivel}|^] T-pyJjJ-Vp) is the line profile 
of the light blocked by the planet centered on v p . Although this 
expression is the result of an expansion in both / and v p , it gives 
a better representation of the measured radial velocity. Hirano 
|et al. ( 2010) 1 also provide more complex expressions in the case 
of Voigt profiles, and in Hirano et al.| ( |201 1) , additional effects 
are taken into account such as macro-turbulence. But the result 
is not expressed as a simple analytical formula and it requires 
several numerical integrations. 

The model of |Hirano etaT1 ( |2010||201 1) > is well adapted to the 
iodine cell technique, but it still does not correspond to the anal- 
ysis routines used on stabilized spectrographs, e.g., SOPHIE and 
HARPS ( TBaranne et aL| [19961 |PeP e et al l |2002"1 >. In these rou- 
tines, the observed spectrum is first cross-correlated with a tem- 
plate spectrum. This gives the so-called cross-correlation func- 
tion (CCF) which can be seen as a weighted average of all the 
spectral lines convolved with a rectangular function. Finally, the 
CCF is fitted by a Gaussian whose mean represents the observed 
radial velocity i>ccf ( see Fig.[T]l. Currently, there is no analytical 
expression of this quantity in the literature. The goal of this pa- 
per is to provide such an expression. As we will see, even in the 
general case where the spectral line profile Tatar is not Gaussian, 
the resulting formula is as simple as Eq. (|2]i derived by Hirano 



et al. (2010 1 and exact in v p . 

This paper is organized as follows. In Sect. [2] we analyti- 
cally derive an unbiased expression i>ccf modeling the RM ef- 
fect as measured by, e.g., SOPHIE and HARPS. This formula 
is specially designed to simulate the radial velocity measure- 
ments obtained by fitting a Gaussian to a CCF. We first provide 
very generic expressions that relies only on the symmetry of the 
spectroscopic lines, and then, we give a much simpler formula 
corresponding to Gaussian subplanet line profiles. In Sect. [3] we 
propose a new expression t^dine of the RM signal derived by the 
iodine technique. In comparison to the previous ones, it is ana- 
lytical and valid for any spectroscopic V sin i A . Then, in Sect. [4] 
we detail the calculation of the parameters entering into our for- 
mulas vqcf and Ui dine> which are the flux fraction / occulted by 
the planet, the subplanet velocity u p , and the dispersion /3 p . In 
Sect. [5] we compare our analytical results with simulated data 
generated with SOAP-T, a modified version of the numerical 
code SOAP (Boisse et al. 2012| l, able to reproduce RM signals 



In this paper, unit Gaussia ns are defined by Q^v) 



exp (-v 2 /(2cr 2 )), while in |Hirano et al.[|2010} , they are de- 
fined by Qe-{v) = ( ^mx) -1 exp (-u 2 There is thus a difference of a 
factor 2 in the parenthesis of Eq. |2]( with respect to ( |Hirano et al.||2010) 
Eq. 36). 
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correlated with 



fitted by 



fitted by 




fCCF 



Fig. 1. Simplified illustration of different methods to compute the Rossiter-McLaughlin effect, vmo, iodine, and uccf represent the 
result of the hypothesis made in|Hirano et al. (2010), the result of the data reduction done with the iodine cell technique, and the 



result of the data reduction done with the CCF technique, respectively. 



(Oshagh et al. 2012[ ). We also analyze biases introduced by the 
application of a wrong model in the fit of RM signals. Finally, 
we conclude in Sect. [6] 

2. Modeling of the RM effect measured by the CCF 
technique 

2.1. General derivation 

We derive a very general expression of the RM effect obtained 
using a Gaussian fit to the CCF. As in Hirano et al. ( 2010| l, we 
only consider the linear effect with respect to the flux ratio /. 
But the method can easily be generalized to higher orders. We 
define fsizriv), Tpi a (v - v p ), and 'Ftransit(f) the line profile of the 
CCF produced by the integrated stellar surface, by the part of 
the stellar surface covered by the planet, and by the uncovered 
stellar surface, respectively. At this stage, the line profiles T sXm , 
Tpia, and transit are not necessarily Gaussian. The dependency 
of Tp\ a is on (u - v p ) because this line is centered on v p . By con- 
vention, Tstuiv) and Tp\ a (v - v p ) are normalized to one. With this 
convention, even absorption spectral lines are positive. We have 
^transit(y) = T^tarO) - /TpiaO - v p )- Moreover, we denote Q^v) 
as the unit Gaussian profile with dispersion cr and centered on 
the origin 



@<r(v) = 



1 



(3) 



The fit of the CCF measured during transit, ^nmsitC^X by a 
Gaussian corresponds to the maximization of the likelihood 



X 2 (d, V, Cr)= J (!F, rans it(» - aQcr(y - vjf dv 



(4) 







da 


-J 




2a 


dv 


~~ ~o^ 


dx 2 


2a 


do- ' 





with respect to the normalization factor a, the mean velocity v, 
and the dispersion cr. The partial derivatives read as 



J(v - v)@a-{v - v)(T lm ns,n(v ) - aQ a {v - vj) dv ; (5) 

J l(v - v) 2 -o- 2 )g a (v-v)(Ttmn S h(v) - a@ a-iv -5)) dv . 

We now set all the derivatives to zero, express Ttmnsni") as a 
function of !F sta i(w) and T-px^v), and reorder the terms. The system 
of equations |5]) is equivalent to 

J F(a, v, cr;v)dv = f J T p ^(v - v p )@ a (v -v)dv; 
J" vF(a, v, cr;v)dv = / J" vT v \Jj> - v p )Q <T {v -v)dv; 
J v 2 F(a, v, cr; v) dv = f J v 2 T v \ d (v - v p )Q a (v - v 



(6) 



) dv . 



with 



F(a, v, cr; v) = Q^v - v)(<F s tm-{v) - aQ^v - vjj . 



(7) 



To solve this system, we apply the usual perturbation 
method. We develop the parameters in series of the flux ratio 
/, i.e., a — aa + fa\ + f 2 a2 + . . ., and idem for v and cr. At the 
zeroth order in /, the system (|6 1 corresponds to the fit of l F s tar(y) 
by a Gaussian. The effect of the planet is absent from the fit, 
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thus vq = 0. The other parameters ao and cr are those of the best 
Gaussian fit outside of transits. 

In a second step, we linearize the system (|6| in the vicinity of 
the zeroth order solution. We thus compute the partial derivatives 
of 



M„(a, v, cr) 



J' 



F(a, v, cr) dv . 



(8) 



appearing in the left-hand side of (Rj), at a — do, v = vo = 0, 
and cr = cr . To avoid fastidious calculation, we make the sim- 
ple hypothesis that Tstta(v) is symmetrical, or more precisely, 
an even function of the velocity v. In practice, due to the con- 
vective blue shift (CB), spectral lines are not perfectly sym- 
metric. Nevertheless, the net effect of CB is to add a constant 
offset that does not modify RM signals (Albrecht et al. 2012| l. 



Then, Mo(a, v, cr) and M 2 {a, v, cr) are even functions of v, while 
Mi (a, v, cr) is odd. As a result, the derivatives dMo/du, dM^jdv, 
dMi/da, and dM\jdcr taken at v — vanish. Moreover, if the 
subplanet line profile is also an even function then the integrals 
of the righthand side of the system (|6]l become simple convolu- 
tions at v = 0. Then there is only 



dM \ (dM \ r _ _ i, , 

= [(Pffar,)*^] (v,,) . 



<dM 



v dv 
'3M 2 



(9) 



)M 2 \ (dM 2 \ v, j„ , _ i, s 

da~ J ° l + \~dcr~ J 0-1 = i ( ^' ro) * Fpla J ' 



where * denotes the convolution product. The first and the third 
lines are independent of the mean velocity v\. Their resolution 
provides a correction to the amplitude and the width of the best 
Gaussian fit during a transit. The second equation is the most in- 
teresting, since it contains the quantity v\ we are looking for. By 
chance, this is also the simplest. The velocity anomaly obtained 
by the Gaussian fit uccf = v = fv\ is then 



y CCF - / 



( dMi 

\~dv~ 



[(f&ro) * !Tpta)] (v P ) 



(10) 



where dM\jdv should be computed at (an, 0, <xo), i.e.. 



3M X 



— - — ;@cr <) (v)( < F s m I (v)-2ao0 cro (v)) dv 
a 



(11) 



4<t yfjr 



+ 



[(^)*(f.t»-flo^)](0) 



The convolution product taken at on the righthand side of (jTTJ 
is nothing else but M%, which cancels by the definitions of ciq and 
ctq. Thus the RM effect now reads as 



VCCF = - 4 °"° — fhvffa-o) * ?>la)l (Vp) ■ 
do 



(12) 



This expression does not depend directly on the stellar line pro- 
file :F star . The dependence only occurs through the best Gaussian 
fit («o and cr ), which is performed to derive the radial velocity. 
The formula < fl2] > is thus very powerful, since it does not require 
any knowledge on Tstsf Unfortunately, the amplitude oq of the 



best Gaussian fit in ( 12 1 is associated to a normalized line profile 
Tstai while, in practice, the area of a CCF is difficult to measure, 



and Tstm is never normalized. We thus provide the expression of 
«o as a function of Tstar and the best Gaussian fit Q^, 



a 



J f?cro(l#r s tar(y) dv 



2cr A^[^ <ro *r star ](0) 



(13) 



We emphasize that ciq is independent of v p . As a consequence, it 
does not affect the shape of the RM effect, but only slightly the 
amplitude (ciq remains close to one). The computation of ciq is 
detailed in Appendix [B| 

2.2. Gaussian subplanet line profile 



At this stage, the expression (|T2J is very general, and holds as 
long as the line profiles T^Av) and T v u{v) are symmetric. 

We now make the hypothesis that the subplanet line profile 
■Fpia(y) is Gaussian. It should be stressed that, as long as the 
planet radius is small compared to that of the star, the subplanet 
line profile is only weakly affected by the stellar rotation (see 
Sect. Q and thus, it is well approximated by that of the non- 
rotating star, which we assume to be Gaussian. Then, the sub- 
planet profile can be considered Gaussian, or equal to the sum 
of two Gaussians if macro-turbulence is taken into account (see 
Appendix [A}, We denote fi p as the width of the subplanet line 
profile, i.e!7jFpi a (i>) = Qp p (v). In that case, the expression of the 
RM effect (fT2ll becomes 



1 

^CCF 

fl 



2oi 



3/2 



fv p exp 



(14) 



This is the main equation of this paper. It represents a good com- 
promise between simplicity and accuracy for the modeling of 
RM signals measured by a Gaussian fit of the CCF. 

2.3. Gaussian stellar line profile 

For completeness, we give the expression in the case where the 
stellar line profile is a normalized Gaussian with dispersion jS*. 
The best fit should give aq = 1 and <t = yS*. Then, we get 



VCCF - - 



3/2 



fv p exp 



m+P 2 P )) 



(15) 



This formula is equivalent to Ohio Q given by Hirano et al. 
( 20 1 0| > . More precisely, Eq. |2]i is the beginning of the expan- 
sion of i>ccf- Indeed, if the stellar line profile is Gaussian, the 
two approaches are identical. 

3. Modeling of the RM effect measured by the 
iodine cell technique 

In this section, we first explain the equivalence between the io- 
dine cell technique and the maximization of the cross-correlation 
C(v) ( |Hirano et al. |2010 2011 1. The aim is to emphasize the 
hypotheses behind this equivalence and, thus, to show its limi- 
tations. In a second step, we provide a general expression that 
models the iodine cell technique. It should be stressed that the 
function C(u) is different from the CCF of the previous section 
and is not used in the same way. It involves the spectrum dur- 
ing transit and a modeled one without transit deformations. This 
function C(u) is computed to provide the RV at its maximum. 
On the other hand, the CCF is the cross-correlation between the 
spectrum and a mask. The goal is to provide a single averaged 
line that is then fitted by a Gaussian curve. 
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3.1. Link between the iodine cell technique and the 
maximization ofC(v) 

The analysis routine based on the iodine cell technique involves 
a fit with 13 parameters of the observed spectrum by a mod- 
eled one that is Doppler-shifted ( |Butler et al.||1996) . We assume 
that this can be approximated by the fit of a single parameter 
(v) representing the Doppler shift between a modeled line pro- 
file Tsiariv) and the observed one (here during a transit) T^ansitC"). 
The chi-square of this fit reads 



X+oo 
(iTtransitO) - T staI (v - 5)) dv . 
OO 



(16) 



The minimization of this chi-square corresponds to dx 2 /dv = 
with 



dx\v) 
dv 



/dT 
-^(v - v)(r aa n S n(v) - T st -Av - v)) dv. (17) 



The integral on the righthand side can be split into the sum of 
two integrals 



dv 



= 2(h{v) + h(v)), 



(18) 



where the first one, I\(~v), is the cross-correlation of dTsm/dv by 
^transit taken at v, while the other is, after the change of variable 
u — v — V, 



h 



-f 



dTsu 
dv 



-(.uWm-Au) du 



(19) 



If y star is even, its derivative is odd, and thus the integral I2 over 
R vanishes. In that case, only 



dx 2 (v) = JdTjto 
dv \ dv 



(20) 



remains, where ★ denotes the cross-correlation product defined 
by 



[f*g](x) 



-f 

•J — CX 



f(y - x)g(y) dy 



(21) 



We then use the property of the derivative of the cross- 
correlation of two functions 



-rif *g) = -~r *g = f * t 

dx dx dx 

By identification, we obtain 



d x 2 (v) 



1 d trzr rjr \ ~dC(v) 

— — 2— (y s tar * /transit) = — 2 — 

dv dv dv 



(22) 



(23) 



where C(v) is defined as in |Hirano et al.lpOlO) . Thus, the min- 
imization of the chi-square involved in the iodine cell technique 
is indeed equivalent to the maximization of the cross-correlation 
C(v) as computed in |Hirano et al. (2010 1. This result holds as 
long as the complicated fit with 13 parameters can be modeled 
by the fit of the single parameter v, and if the modeled line profile 
is symmetrical. The second condition may not be true in general. 
If the asymmetry is not too strong, the integral I2 would be a 
small perturbation, and the result obtained by the maximization 
of the cross-correlation C(v) should differ from the minimization 
of the chi-square by only a small constant. 



3.2. General expression of the RM effect 

To derive a general expression of the RV signal measured by the 
iodine technique, we use the same model as Hirano et al. ( 2010) >, 
which consists in maximizing C(v) where 

Civ) = J T sVdI (v - v)[r stax (v) - fTfrQ) - v p )] dv . (24) 

Then, the condition dC(v)/dv = leads to 



d^sta 

dv 



*r sti 



® = / 



dT sti 
dv 



-(v-v)T ? i d (v-v p )dv . 



(25) 



As in the previous section, we expand v in series of /: v = vq + 
/Di + .... At the zeroth order, we get 



dv 



Cvo) - 



(26) 



If ? r star(y) is an even function, this equality gives Do = 0. 
Otherwise, Do would be a small constant, depending only on the 
shape of the spectral lines, but not on v p . Here, we assume that 
Do = 0. At the first order in the flux ratio /, assuming that Tpia is 
even, we obtain, with iodine = fvi , 



- If 

^iodine ~ a J 



dv 



phi 



(Vp) 



where 



d 
dv 



dT sta 
dv 



dv 



dv. 



(27) 



(28) 



This result is more complex than ( |T2] i because it involves the 
derivatives of the stellar line profile Tstar instead of the deriva- 
tives of a best Gaussian fit Q aa which are analytical. Of course, 
if ^tar is Gaussian, we retrieve the result of the previous section 
(see Sect. |2.3[ >. 

In Appendix |C| we provide an analytical expression of the 
numerator of t^dine ( 27 1 in the case where the subplanet line pro- 
file is Gaussian and for a general limb-darkening law. 

4. Parameters of the subplanet line profile 

The expressions of the RM effect ( fl4| and ( |27") > depend on the 
fraction / of the flux covered by the planet, the subplanet veloc- 
ity v p , and the dispersion /3 p . There are two approaches to evalu- 
ate them. On the one hand, both the flux fraction / and the mean 
veloci ty v p can b e computed exactly as a series of Jacobi polyno- 
mials (Gimenez 2006 1. This is useful in the case of binary tran- 
sits where the occulting object is big. On the other hand, only 
the flux fraction / is derived exactly using analytical algorithms 
such as the one given by (Mande l~& Agol| |2002[l , while v p and 
Pji are estimated assuming uniform intensity bel< 



intensity below the planet 
( |Hirano et al.| |2010] |201 \) . Then, if (x, y) are the averaged co- 
ordinates over the surface of the star covered by the planet, and 
normalized to the radius of the star, v p w xV sin/*, while j3 p is 
constant and represents the width /3q of the nonrotating star line 
profile. 

Here, we choose a compromise between the two approaches 
and take the slope and the curvature of the intensity below the 
planet into account. This gives a better estimate of v p and fi p in 
comparison to the uniform subplanet intensity hypothesis. But 
also it turns out that the method provides a simple and accu- 
rate expression for the flux fraction /. Another advantage of this 
method is that it can be easily applied to more complex problems 
where the gravity-darkening or the tidal deformations of both the 
planet and the star are taken into account. 
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4.1. Method 



4.2. Subplanet velocity 



To describe the method, we take the example of the computation 
of the flux fraction /. The expression of / reads as 



I I I(x, y) dx dy , 
JJs„ 



The above method applied to I(x,y) to get the flux fraction / 
can be adapted to any other function. For example, the subplanet 
velocity is defined by 



(29) 



v p = Vsin/* 



where S p is the surface of the stellar disk covered by the planet 
normalized by the square of the radius of the star. I(x, y) is the 
normalized limb-darkening of the star expressed as a function 
of the normalized coordinates (x, y). For the moment, we do not 
need to give its expression. 

A very rough approximation of / is obtained assuming uni- 
form intensity below the planet. In that case, we obtain 



ff s xl(x,y)dxdy 
ff s I(x,y) dxdy 



(36) 



We denote H^, H^J, and as the components of the Hessian 
of xl(x, y) at the averaged position (x, y) ( [3T) . Using the expres- 
sion of / ( 35 1, at first order in r 2 , we get 



/ = I(x,y) jj^ dxdy 

= I(x,y)S p , 
where 



V sin L 



= x + 



(fl2>-*H£)(C*-tf) 



(30) 



2I(x,y) 

+2(H^-xH^)((x-x)(y-y)) 



(37) 



= — — \\ xdxdu , and u = — — \\ udxdu , 



(31) The denominator I(x, y) in d37b, as well as the terms in H\ ._, . , H, 



and , come from the expansion of the denominator of (|36|). 



j(0) iA0) 
iiiy > 



are the coordinates of the barycenter of the portion of the stellar 
disk covered by the planet. The formula d30l) works well for very 



small planets but only during full transits (e.g. |Mandel & Ago! 
2002). Close to the limb, the intensity I(x, y) varies strongly with 



position and this approximation is not valid anymore. To over- 
come this issue, we propose to make an expansion of the limb- 
darkening profile I(x, y) in the vicinity of (x, y). We get 



4.3. Width of the subplanet line profile 

The width of the subplanet line profile f3 p is a combination of the 
width of the nonrotating line profile /3o and a correction 6/3 p due 
to the rotational broadening 



J- = /(*, y) <1> + Jf> (x-x) + Jf (y-y) + ±H™ ((a 

+±H® ((y - yf) + H% ((x - x)(y - y)) , 
where for any function f(x, y), 

(f(x,y)) = 4- \\ f{x, y) dxdy, 



P 2 p =/3 2 o + W 2 p 



(38) 



We define 'Vi as the average of the square of the subplanet ve- 



(32)locity 



T 2 = (VsimV) 2 



(33) 



where jf^ = d x I(x, y) and jf^ = d y I(x, y) are the components 
of the Jacobian of the surface intensity I(x, y) computed at (x, y). 
Similarly, the components of the Hessian are 



$ s x 2 I{x,y)dxdy 
jj s I(x, y) dx dy 

With this notation, we have 
6pL =<V 2 - vl . 



(39) 



(40) 



H' 



(0) 



8 2 I(x,y) 



H (Q) 

yy 



d 2 I(x,y) 



H (P) 

xy 



d 2 I(x,y) 



Noting H^, Hfy, and H®, the components of the Hessian of 
x 2 I(x, y) at the averaged coordinates (x, y), the expression of ^2 
reads as 



(34) 



dx 2 m dy 2 dxdy 

By construction, x and y are defined by x — (x), and y = (y). 
Thus, the linear terms in factor of the Jacobian (J x °\ jf^) cancel 
in (32 1. At that point, only 



<V 2 



(Vsin/*) 2 



= i 2 + 



1 



(Hg-#H%)(<?-tf) 



2I(x,y) 



(41) 



/ 



I ix ,y )+ ^((x-x) 2 ) + ^((y-y) 2 ) 



+2(fl« 



2 H^)((x-x)(y-y)) 



(35) 



+H^((x-x)(y-y)) 



remains. The first term i n ([35} corresponds to the rough approx- 
imation derived in Eq. ( 30 1. The other terms provide a correc- 
tion proportional to the square of the normalized planet radius 
(r = R p /R+) and are expected to be small. 



At first order, the i 2 in 'Vi (39 1 cancels with the square of x in 
the expression of v p (37 1. Thus, at first order, S/3 p vanishes and 
the width of the subplanet profile is equal to the width of the non- 
rotating starjS ( ). However, the quadratic terms do not cancel, and 
this provides an estimation of the contribution of the rotational 
broadening to the actual width of the subplanet profile. 
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Table 1. Second derivatives of x"I a (x, y), Eq. (42 1. 



H%> = n(n - l)x"- 2 4(x, y) - n(n + l)ax"T a . 2 (x, y) 
+a(a - 2)x" +2 I a - 4 (x, y) 

= -nax"I a - 2 (x, y) + a(a - 2)x"y 2 I a - 4 (x, y) 



h: 



—nax" yl a -2(x, y) + ct(a — 2)x l,+> y l a - 4 (x, y) 



4.4. Limb-darkening and its derivatives 

As we saw above, the determination of the subplanet profile 
depends on the limb-darkening law and its second derivatives. 
In this section, we provide generic formulas assuming that the 
limb-darkening law is a linear combination of functions I a (x, y) 
defined by 



(42) 



y is the cosine of the angle between the 



I a (x,y)=fi a = (l -x : 

where /j = -y/l - x 2 - 
normal of the stellar surface at (x, y) and the observer. 

The second derivatives H^, H y " y \ and of x"I a (x, y) are 
given Table[T] These are the ones that are needed to compute the 
flux fraction /, the subplanet velocity v p and the dispersion f} p . 
In practice, only the cases n — 0, 1, 2 are used. 

From these general formulas, one can derive the expressions 
for the quadratic limb-darkening which reads as 

/,(*, y) = I q (0) (l - m(l - M ) - " 2 (1 - V?) 

with I q (0) the central intensity such that I q (x, y) is normalized to 
one. By identification, we get 



1 — U\ — U2 



11 7T(1 - Ml/3 - M 2 /6) 

u\ + 2t<2 



M, = 



n(l — Mj/3 — u.2/6) 

~U2 



(44) 



n{\ — «i/3 - M2/6) 
Equivalently, the so-called nonlinear limb-darkening is usually 
expressed as 



I n i(x, y) = l„i(0) 



i~2>(i-^ 2 ) 



with 



C n = 



1 - Ci - Ci - C3 - C4 



tt(1 -ci/5-c 2 /3-3c 3 /7-c 4 /2) ' 
and for 1 < n < 4, 



(45) 



(46) 



(47) 



" tt(1 -CJ/5-C2/3-3C3/7-C4/2) 

The normalizations in ( |44] >, |46|, and ( |47| i have been deduced 
from the integral of each I a (x, y) over the entire disk of the star 

(48) 



J J //" dxdy — J dcf> J 



rf0 cos ff+1 6»sin6» : 



a + 2 




Fig. 2. Definition of the angles tp p l \ ip p 2 \ ip^, and ip^ during 
partial transit. The large circle centered on S represents the star, 
and the smaller one, centered on P, is the planet. 



4.5. Averaged coordinates and covariances 

The last quantities that need to be computed in order to get 
the subplanet line profiles are the averaged coordinates (x), 
(y), the variances ((x - x) 2 ^, {iy-y) 2 ^, and the covariance 
((x - x)(y - y)). For that, we distinguish two cases. 



4.5.1. During a full transit 

In the case of a full transit, i.e., when the disk of the planet is 
fully inside the disk of the star, the problem gets simpler since 
the integrals (33 i have to be computed over a uniform disk of 
area S p = nr z and centered on the coordinates (xq, yo) of the 
planet. We get (x, y) = (xq, yo), and 



(49) 



((x-x) 2 )=((y-y) 2 ) = r - , 
((x - x)(y - y)) = . 

4.5.2. During ingress or egress 

If the disk of the planet is crossing the limb of the star, the area 
where the integrals of the form ( 33 1 are computed is not circular 
(see the shaded area i n Fig.[2| . In that case, we use the very pow- 
erful method of |Pal 
work also for mutua 



(2012), which gives expressions that also 
transits. 

We recall briefly the method that relies on Green's theorem 
converting an integral over a surface into an integral over the 
contour of that surface: 

oj(x, y) dxdy = (j) (f x {x, y) dx + f y (x, y) dyj . (50) 

In this equation, dS is the boundary of S, and uix, y) the exterior 
derivative of f(x,y) = (f x ,f y ) defined by 

dfy 8f x 



co(x, y) 



dx dy 



(51) 
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When the planet is crossing the limb of the star, the boundary is 
the union of two circular arcs. One of them follows the edge of 
the planet centered on (xo, yo) with radius r. The coordinates of 
any points of this arc and the tangent vectors are of the form 



(52) 



x — xq + r cos <p , dx — —r sin tp dtp , 

y = yo + r sin ip , dy — r cos tp dip . 

The angle ip varies between two limits tpy and tpy , correspond- 
ing to the intersections A and B between the circumferences of 
the planet and of the star, respectively (see Fig. [2). The second 
arc fellows the edge of the star and is parameterized by the co- 
ordinates 

x = cos w , dx — — sin <p dip , 

Y * * (53) 

y — sin ip , dy — cos ip dip , 

with tp going from ip\^ to tp^ associated to the intersections B 
and A, respectively. 

More generally, if we denote (xj, yi)j=p, s as the centers of the 
arcs, and (rj)j =ps as their radii such that 

and y(ip) - yj + rj sin <p , (54) 



x((p) — xj + rj cos tp , 
we obtain ( [PSl] [2012] ) 

u>(x, y) dx dy — 



ffs,, 

(55) 

(fy(x(<p), y(<p)) cos <p - f x (x(tp), y(ip)) sin ip) rj dip . 

j-p* -J 

Here, we are interesting in the cases where u>(x, y) stands for 
1, x, y, x 2 , y 2 , or xy. The field vectors f(x,y) associated to 
those a>(x, y) are not uniquely determined. We choose (-\y, \x), 
(—xy, 0), (0, xy), (-x 2 y,0), (0,xy 2 ), and (-\xy 2 , \x 2 y), respec- 
tively. We denote o>y = x l y' and f f j(x, y) as the functions whose 
exterior derivative is o>y. The integral of f f j along a circular arc 
reads as 



Fiji?) 



f(f„ 



,(xq + r cos ip, yo + r sin <p) cos <p 



(56) 



-fij, x (xo + r cos tp, j/o + r sin tp) sin tp^rdtp . 

Since the /y are polynomials in x and y, the F,-,- c an be com- 
puted using the recurrence relations provided by Pal (2012). The 
results are displayed in Table [2] for i = 0,1,2, and j = 0,1,2. 
From them, one can derive the quantities present in the expres- 
sions of the flux fraction /, the subplanet velocity v p and the 
width/?,,: 



^01 



(57) 



O)20 - X 



o)oz - y 



((x - x)(y - y)) 



Snff 



xy . 



1 1 1 

Mandel & Agol ° 
ui = 1, U2 = — 
u\ = 0, U2 = 1 — 




O 

X 

< 

-10 

.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 
phase 

Fig. 3. Transit light curves for r = 0.1 and quadratic limb- 
darkening with two sets of coefficients u\, U2- The solid curves 
in red and in green are obtained from the approximation ( [3~5| , 
while the black open circles are computed using the routine of 



Mandel & Agol (2002 1 



5. Comparison with simulations 

5.1. Transit light curve 

Although it was not the main goal of this present work, in the 
derivation of a precise modeling of the RM effect, we obtained a 
new expression of the flux fraction / occulted by a planet during 
a transit (see Eq. (35 i). In comparison to existing formulas that 

, the one of this 



are exact (e.g. Mandel & Agol 



2002 Pal 



2012) , 



paper relies on an expansion of the intensity in the vicinity of the 
averaged position of the planet. We thus expect our formulation 
to be less precise. 

Figure [3] shows the comparison between the approximation 
( |35) and the exact formula derived by Mandel & Agol ( 2002). By 
eye, it is not possible to distinguish between the two approaches. 
In the residuals, however, we can see that the maximum of de- 
viation occurs close to the limb, more exactly, when the edge of 
the planet is tangent to that of the star. Indeed, at the border of 
the star, the limb-darkening becomes steeper and steeper, and the 
derivatives d x I a (x, y) and d y I a (x, y) even go to infinity for a < 2. 
Nevertheless, this singularity is smoothed out by the decrease in 
the overlapping area between the planet and the star disks during 
ingress and egress. 

One advantage of the present formula is that it can be easily 
generalized to more complex problems, as in the cases of a dis- 
torted planet, distorted star, important gravity limb-darkening, 
and so on. For our purpose, it provides an accurate enough esti- 
mation of the flux that can be used to derive the RM effect. 



5.2. Subplanet profile 

We checked the accuracy of our new formulas of the subplanet 



velocity 



Q and the width B 2 p = + SB 2 with 5B p given 
by (BO). For tnat, we used the software called SOAP, for Spot 
Oscillation And Planet (Bois se et al.||2012) , to produce artifi- 
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Table 2. Integrals Fy (56 1 of the field vectors /y(jc, y) along a circular arc centered on (xq, yo) with radius r. 



Fo,o(v) = x + x ° sin f ~~ y cos 

FijoQp) = ~ x oyo rcos *P - ^fo'' 2 cos2 f + ~^ x oi" 2 (2<p - sin 2<p) + y^ r3 (3 siny? - sin 3ip) 

1 , . 2 1 2 . 1 3 

Fo,i(</>) = xo!/o''sin</? + -xor" sin <p + -^yor (2ip + sin2</>) — — r (3 cosy + cos 3ip) 

Fioif) - -x^yor cos ip - x y r 2 cos 2 <p H — xlr 2 (2ip - sin2</j) yor 3 (3 cos ip + cos 3<p) H — x r 3 (3 since - sin 3cc) H r 4 (4cp - sin4cp) 

4 12 6 32 

pQiif) - ■Wo'" since + -fof/o'' 2 sin 2 <P "t — yl^i^-V + sin2co) H x r 3 (3 since - sin3cp) j/or 3 (3 cos c/> + cos 3<p) H r 4 (4co - sin4cp) 

4 12 6 32 

Fi,i(<P) - - x oyo r (2r<p + x o sin<p - y Q cosy) -\ — x\r 2 sin 2 ip (y 2 , + r 2 )^ cos 2 <p H t/o^C 1 5 siny - sin3y) x r 3 (l5 cos ip + cos 3<p) 

4 8 8 48 48 



cial data as close as possible to real observations. This code is 
a numerical tool that models radial velocity and photometry ob- 
servations of stars with spots. It has been updated recently to 
also model the effect of a planet transiting a spotted star, and 
was renamed SOAP-T (Osh agh et al.| 2012| l. Briefly, the code 
divides the disk of the star into a grid. To each cell of that grid, 
a Gaussian profile with a width /?o and amplitude I(x, y) (in our 
notation) is assigned. This represents the intrinsic line profile of 
the nonrotating star as detected by the instrument. These lines 
are then shifted in velocity according to their position with re- 
spect to the spin-axis and the V sin i+ of the star. All the lines of 
the cells uncovered by any spots or planets are added together 
to produce an artificial CCF that is then fitted by a Gaussian to 
derive a radial velocity. 

With SOAP-T, we produced the CCF of a star with a tran- 
siting planet at different positions of the planet on the disk. We 
also generated the CCF of the same star while the planet is not 
transiting, and by taking the difference, we got the subplanet pro- 
file. Such profiles are displayed in Fig. [4]for different values of 
V sin i+. Unless specified explicitly, here, and in all the following 
simulations, the star is a solar-type star with a quadratic limb- 
darkening law whose coefficients are u\ = 0.38, 112 = 0.3, and an 
intrinsic line width without rotation of ySo = 3 km.s -1 . The planet 
is a Jupiter evolving in the equatorial plane of its star, its radius 



is r = Rp/Rseu = 0.1099. In Fig. 
low rotating stars are Gaussian. 



41 the subplanet line profiles of 
s results from the hypothesis 



of SOAP-T, which assumes Gaussian intrinsic line profiles. But 
we observe that the Gaussian shape holds even for V sin/* = 20 
km.s -1 , which validates our assumption leading to Eq. ( 14 1. 



To each of the artificial subplanet profiles generated with 
SOAP-T, we also computed the mean velocity v p and the disper- 
sion f3p, to be compared with our formulas ( 37 1 and ( |40| >. Figure 
|5]shows the results for v p after normalization to remove the effect 
of the V sin i+ of the star. We checked that the figure is indeed un- 
changed up to V sin/* = 20 km.s -1 . The numerical outputs ob- 
tained with SOAP-T are plotted against two different analytical 
approximations denoted So and S2. In Su, the surface brightness 
of the star is taken uniform below the disk of the planet, while 
in S 2, the second derivatives are taken into account as in (37) . 
We observe that where the error is maximal, close to the limb, 
S 2 improves the determination of v p by about a factor 3 with re- 
spect to So- In the case r - 0.1 and V sin/* = 10 km.s -1 , the 



O 



CD 

3 



o 

r—i 

X 



SOAP o 

20 km/s ■■■ 

10 km/s — 

1 km/s — 




velocity (km/s) 

Fig. 4. Example of subplanet line profiles obtained with SOAP- 
T (circles), compared with Gaussian profiles (curves) for differ- 
ent stellar V sin/*. 



maximal error provided by S 2 is about 20 m/s which represents 
a relative difference of 0.2%. 

In the case of the dispersion {} p , the difference between the 
estimation derived assuming uniform (So) and nonuniform (S 2) 
brightness below the planet disk is more evident (see Fig. 16). 
Indeed, in the former case, /3 P remains constant and equal to the 
width Po of the nonrotating star line profile, while we observe 
that for the simulated and the modeled line profiles, the shape 
of f3 p as a function of the orbital phase looks like a trapezoid 
with the large base at ySo and the maximum at approximately 



4 



'y82 + ( r ysin/*) 2 /4. 
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Fig. 5. Subplanet velocity v p produced with SOAP-T (blue 
points), approximation So assuming uniform intensity below the 
disk of the planet (red curve), and approximation S2 taking the 
second derivatives of the stellar surface brightness into account, 
Eq. ( 37 1 (green curve). 
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Fig. 6. Subplanet dispersion f) p produced with SOAP-T (blue 
points), approximation So (red curve), and approximation S2, 



Eq. (40 1, (green curve). 



5.3. Rossiter-McLaughlin effect 

We now compare our analytical expression of the Rossiter- 
McLaughlin effect uccf ( 14 1 with signals generated with SOAP- 
T, which simulates the reduction analysis of the CCF technique 
numerically. 

Figure [7] displays the results for different Vsin/*. As long 
as Vsin/* is below or equal to 10 km.s -1 , the error induced 
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Fig. 7. RM signals produced with SOAP-T (solid black curves) 
for different Vsin/*, and results of uqcf ( p"4"l ) (different dashed 
curves). 



by the analytical formula remains lower than ~1 m/s, which is 
close to the magnitude of the precision of RV measurements. 
In that case, the analytical approximations are almost indistin- 
guishable from the numerical simulations. However, for larger 
V sin the agreement between numerical signals and analytical 
ones is weaker. For example, when V sin i* = 20 km.s -1 , the an- 
alytical approximation leads to a maximal error of 10 m/s, which 
is 5% of the amplitude of the signal. Nevertheless, it should also 
be noted that for fast-rotating stars the spreading of the spectral 
lines over the detectors decreases the precision of the measure- 
ments. In any case, the analytical expression uqcf brings a defi- 
nite improvement over other formulas, which have not been de- 
signed to simulate the CCF technique as we see in the following 
section. 

5.4. Comparison between different techniques 

To highlight the effect of the instrument and of the data re- 
duction analysis, we generated different models of line profile 
and compared the RM signals computed numerically with the 
results of the analytical formulas dccf (14 1 and undine (27 1. In 
our examples, the line profiles are of three types: :Fhires(c), 
y r HARPs(f), and :Fcoralie(u)- They are associated to three RM 
signals: uhires, Sharps, and Ucoralie, respectively. It should be 
stressed that the goal is not to reproduce the lines observed by 
those instruments exactly, but to capture their main characteris- 
tics. On the one hand, HIRES and HARPS are two spectrographs 
with high resolutions that we assume to be identical with a width 
fio = 2.6 km.s 1 for nonrotating solar-type stars. On the other 
hand, the resolution of CORALIE is about twice lower, and the 
intrinsic w idth of the same stars is about fio = 4.5 km.s -1 (Santos 
etal. 2002b. We consider a star with V sin/* = 15 km. s~', which 



is adapted to our illustration. Finally, the transiting planet is a 
Jupiter-like planet with a radius R p = 0.1^ star . 

Figure [8] shows the simulated line profiles. The panel [8^ 
displays the rotation kernel Hiv), a stellar profile Tm m {v) with 
/?o = 2.6 km.s -1 , and a subplanet profile T v \Jj) ~ Vp) multi- 
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modeling an average line profile observed with HIRES in solid black, and a CCF observed with HARPS in dotted cyan. The same 
with j6o = 4.5 km.s 1 and /3 p - 4.56 km.s 1 represents a CCF observed by CORALIE, in dashed violet. 



plied by the flux fraction /, and computed with f} p — (0^ + 
(rV sin ; + /4) 2 )^ 2 . Figure [8(3 depicts the resulting line profiles 
during transit Twkes(v), Sharps (^)? and !Fcoralie(^)- Following 
our hypothesis, ^FhiresW is identical to Sharps 

From the simulated line profiles, we derived RM signals nu- 
merically. The signal i>hires was obtained from Thires using 
the iodine cell technique, i.e., by fitting the best Doppler shift 
between a line without transit deformation, and the line profiles 
computed during transit. Both Sharps, and i>coralie are the re- 
sults of applying the CCF technique, i.e., a numerical fit between 
a shifted Gaussian and Sharps and ^coralie, respectively. We 
also generated i^max by maximizing the cross-correlation C(v) 
between the line profiles Shires at and out of transit. These four 
RM signals are represented in Fig. [9^. It is notable that the RM 
effects associated to the three instruments are all different. The 
variation between Sharps and ucoralie is on ly due to the change 
in resolution. However, in the case of i>hires and Sharps, the 
simulated lines are exactly identical. The observed difference in 
the RM signal is the result of the chosen data reduction tech- 
nique. Figure [9^ also confirms that the maximum of the cross- 
correlation C(D) gives the same result as the iodine cell technique 
(when the stellar lines are symmetrical) since i>hires = ^cmax- 

The last three panels of Fig. |9]represent the comparison be- 
tween the simulated RM signals and the analytical formulas uccf 
( |14) , and iodine p7) associated to the CCF and the iodine cell 
technique, respectively. We observe that the formulas adapted 
to the analysis routines are in good agreement with the respec- 
tive simulations. We also notice that for CORALIE, whose res- 
olution is lower, the two analytical formulas give roughly the 
same result. This is because the stellar line is less affected by 
the rotational kernel and is more Gaussian. We show that, in that 
case, the two methods should indeed provide the same result (see 
Sect.|23). 

From this study, we conclude that a given star observed by 
two different techniques should present two distinct RM signals. 
To date, this notable result has not been seen since the instru- 
ments with the highest signal-to-noise, HIRES and HARPS, are 
located in two different hemispheres. This makes it difficult to 
observe the same stars. For those observed with other instru- 



ments, the expected gaps are diluted by the measurement un- 
certainties. Nevertheless, with the advent of HARPS-North, we 
may observe such discrepancies in the future. 



5.5. Biases on fitted parameters 

As a final test, we simulated artificial data from either the CCF 
or the iodine cell model, and we fit each of these data with the 
two models separately. The goal is not to perform an exhaus- 
tive study of the biases introduced by the application of a wrong 
model in the process of fitting data, but to give an example with 
some typical parameters. 

For this illustration, we considered only one set of parame- 
ters. As in the previous section, the star has V sin i+ = 15 km.s -1 
with intrinsic line width y6o = 2.6 km.s -1 . We chose a quadratic 
limb-darkening characterized by u\ — 0.69 and U2 = 0.0. The 
planet's radius is taken equal to R p = 0AR stal . The impact pa- 
rameter of the orbit is assumed to be 0.3/? star . For information, 
this value is that of a planet with a semi-major axis a = 4/? star 
and an inclination i = 85.7 deg. All these parameters were fixed 
throughout all the simulations. Only the projected spin-orbit an- 
gle /li n p U t was varied from to 90 degrees by steps of ten de- 
grees. For each value of ^i npu t, and each model, 1 000 datasets 
were generated with a Gaussian noise of 10 m/s. Each simula- 
tion contains 50 points, among which 32 are inside the transit 
and 18 outside. In each case, we fit both the V sin /* and the pro- 
jected spin-orbit angle. 



Figure[T0]shows the results of this analysis in the case where 

As 



14 1) 



the data are generated with the CCF model uqcf (Eq 
expected, the parameters recovered with the appropriate model 
are accurate, while those deduced from the iodine cell technique 
formulas are biased. The bias on V sin /* is systematically pos- 
itive and also the most important, especially at large projected 
spin-orbit angle (A invul » 90 deg) where we get (Vsin/*) fit = 
20.7 + 0.5 km.s -1 instead of 15 km.s -1 . One can notice that this 
agrees with the results of, e.g., Simpson et al. ( 2010 ) who applied 
the model of lHirano etakl on WASP-3 observed with SOPHIE. 
They fit a Vsin/* = 15.7*} 1 km.s -1 while the spectroscopic 
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Fig. 9. Comparison of a simulated RM signal observed by different techniques and/or instruments. We used the line profiles of 
Fig. [8]}. The open diamonds, circles and squares represent the RM signal obtained numerically with the iodine cell technique on the 
HIRES line profile, and with the Gaussian fit to the HARPS, or CORALIE, CCFs, respectively. In (a) The gray curve corresponds 
to the numerical maximization of the cross-correlation C(u) of the HIRES profiles inside and outside transits. In (b), (c), and (d), the 
numerical RM signal computed on the HIRES, HARPS, and CORALIE line profiles, respectively, are compared with the analytical 
formulas i>ccf ( 14 i in solid red and iodine (27 1 in dash-dotted green. 



value is only 13.4 ± 1.5 km.s -1 . On the other hand, the bias on 
the fitted projected spin-orbit angle Af\ t remains within 2-cr. This 
parameter is thus less affected by the model. 

The difference in behavior between (Vsin; + )fi t and A& t is 
more evident in Fig. [TT] In that case, the data were simulated 
with the formulas associated to the iodine cell technique: iodine 
(Eq. (27 1). As in the previous test, using the same model for both 
the generation of the data and the fit, leads to very accurate es- 
timations of the parameters, while the application of the wrong 
model introduces biases. The (V sin;*)^ is systematically nega- 
tive as we could expect since the situation is the opposite of the 
one in the previous paragraph. Nevertheless, the error is smaller. 
In the worst case, we get (V sin/+)fit = 11.9 + 0.3 km.s -1 , which 
represents on error of the order of 3 km.s 1 , while it was almost 



6 km.s -1 in the previous example. The situation is similar for A^. 
We observe small biases anticorrelated with those of the previ- 
ous test, but now, the inaccuracy remains within 1-cr. 

We stress that we only fit two parameters in this study, while 
the others are fixed to their exact values. We already observe 
that the best fits tend to compensate for inaccurate models by 
introducing biases. With more free parameters, there are more 
possibilities to balance the model, and it is thus difficult to pre- 
dict the behavior of the fit. Since the models are not linear, we 
should expect the presence of several local minima. Eventually, 
in some of them, the projected spin-orbit angle might be more 
biased than in our tests. This should be analyzed on individual 
case bases, which is not the goal of this paper. 
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Fig. 10. Results of the fits with two different models of mock 
data generated by the CCF technique formulas. In each panel, 
the area with the lightest color represents the two-sigma limit, 
the darkest color is the one-sigma threshold, and the curve is the 
best value. 



fort was made to model the effect of the rotation of the star on 
the width of the subplanet line profile. Since our formulas do not 
rely on any expansion in powers of the subplanet velocity u p , our 
results remain stable even for fast-rotating stars. 

The advantage of having a purely analytical expression to 
model the RM effect is the rapidity of computation. It can thus 
be used to analyze a large sample of RM signals uniformly. As 
a complement to this paper, we make our code accessible to 
the community as a free open source software package. This is 
a library called ARoME, an acronym for Analytical Rossiter- 
McLaughlin Effect, designed to generate analytical RM signals 
based on the formulas of the paper. It also includes the effect of 
macro-turbulence as described in the Appendix |A| The library 
provides a C interface and, optionally, a Fortran 77 or 2003 inter- 
face to be called by an application. The fully documented pack- 
age can be downloaded from the webpage www . astro . up . pt/ 
resources/arome 

Besides the modeling of the RM effect, we also analytically 
derived a new expression for transit light curves (35 1. Although 
this expression is the result of a Taylor expansion of the intensity 
and is only adapted to small planets, it gives good approxima- 
tions. Moreover, the expression is general enough to be easily 
extended to more complex problems. 
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Fig. 11. Same as Fig. 
cell technique formulas. 



6. Conclusion 

One of the main objectives of this paper has been to highlight 
that there is no unique way of measuring RM effects and that dif- 
ferent techniques provide different values of RV anomalies. RM 
signals should thus be analyzed using the appropriate model to 
avoid any biases, at least in V sin/*. This is particularly impor- 
tant in the case of low-impact parameters (planet passing close 
to the center of its star) since then, the projected spin-orbit angle 
only depends on the amplitude of the RM signal. 

We provided a new analytical formula specially designed to 
model RV anomalies obtained by fitting a Gaussian function to 
the CCF, as in the analysis routines of HARPS and SOPHIE. 
We also revisited the modeling of the iodine cell technique, as 
used with HDS and HIRES, for which we derived an analytical 
expression adapted to non-Gaussian stellar line profiles. An ef- 
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Appendix A: Macro-turbulence 

Here, we study the effect of macro-turbulence on the Rossiter- 
McLaughlin signal. We consider only the "radial-tangential" 
model as in ( |Hirano et al.| |201 In that case, if we denote 
Tq(v) as the line profile of the nonrotating star without macro- 
turbulence, the subplanet line profile reads as 



(A.l) 



where M(v) is the rotation-turbulence kernel given by (Gray 
[2005] ), 



M(v) 



and 



JJ" I(x, y)&(x, y, v — xV sin /*) dxdy . 



(A.2) 



&(x, y, v) = i (® R (x, y, v) + & T (x, y, vj) 
- e -(fib) _| 



1 



1 



(A3) 



2 Viu cos 6 £sin#~ 

We highlight the different dependencies on (x,y), on the one 
hand, through cos 6 — y 1 - x 2 - y 2 and sin 9 = ^x 2 + y 2 , and 
on the velocity v, on the other. The coordinates (x, y) are normal- 
ized by the radius of the star. Since is the sum of two Gaussians 
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® R and &t associated to the radial and the tangential broaden- for the tangential broadening. We thus have 
ings, respectively, we also split M(u) into two parts 



tR) 



M(v) = -(M R (v) + M T (vj) , 



(A.4) 



such that M R (u) is associated to ® R , and M T (v) is associated to 

® T . 

Now, we compute the moments of (Mj(u))j =R j as in Sect. [4] 
to evaluate the effect of the rotation-turbulence kernel on the sub- 
planet profile. We have 



(A.5) 



1 f + °° 

A j J-OO 

x J"J~ dxdyl(x, y)®j(x, y, v — xV sin /*) , 
where Aj is a normalization constant whose expression is 

X+OO 
du II dxdyl(x, y)®j(x, y, v - xV sin /+) . (A. 6) 
OO *J*JSp 

Inverting the integrals, we get, for the normalization, 



with 



£ + <V 2 and <V { p = ti+<V 2 

ff s cos 2 8 1(x, y) dx dy 
ff I(x,y)dxdy 

^ff s sin 2 6 1(x, y) dx dy 
ff s I(x,y)dxdy 



(A. 14) 



(A. 15) 



(A. 16) 



where corresponds to the case without macro-turbulence 

c 2 + c 2 



(Eq. \39\). It should be noted that &+& = £ 2 . Let /? be the 



dispersion of To(v), and <5/3 p the dispersion due to the rotational 
broadening alone ( |40| . The subplanet line profile T p i a ( A. 1 1 can 
be approximated byThe sum of two Gaussian functions 



T v \ a (v - u p ) = ]^{Qp R (v - u p ) + Q Pt {u - u p )j 
centered on the same value v p with respective dispersions 

r*r* /~»+oo 

Aj= II dx dy I(x, y) I du ©/*, y, v - xV sin i*) . (A.7) $\ = $ + 6/3p + $ , (A. 1 8) 

•JUSn xj —OO 



(A. 17) 



Since ®j(x, y, u) (A.3i is normalized, the inner integral on the 
velocity is one. It thus remains only 



and 

/3 2 T = f3 2 ) +6f3 2 p +£. 



(A. 19) 



A R =A T =A 



I(x, y) dx dy . 



(A 8) ^ n tn * s ex P ress i° ns ' Sfipi £k> an d Ct are functions of the position 
of the planet on the stellar disk. With this model, the Rossiter- 



McLaughlin effect, as measured by the Gaussian fit of the CCF, 

as in Sect. B] We now focus on the numerator of (A.5 1. By con- re ads as 



struction, y) M - L Then, using the inversion of integrals, we 



get 



uccf - - 



dy 



2fln 



I(x, y) 



X 



X+oo 
du u ®j(x, y,u - xV sin ;*) 
OO 



(A.9) 



( 2cr 2 } 


3/2 

ftp 


exp 


f 


v 2 \ 

p 








( 2 -o 1 


3/2 

fVp 


exp 


/ 


»l } 


W 2 Q +/3 2 T ) 


\ 


2(cr 2 +/3 2 T )j 



(A.20) 



The inner integral over the velocity u gives xVsin/*, we have 
thus 



Appendix B: Normalization factor of the Gaussian 
fit 



Vp = Vsim* 



xl(x,y)dxdy 
ff s I(x,y) dxdy 



for each broadening: radial (j = R) and tangential (J = T). This 
is identical to ( 36 1. Finally, the second moment reads as 



B. 1 . Without macro-turbulence 

(A. 10) Here we detail the computation of the amplitude a of the best 
Gaussian fit ( 13 1. In a first step we neglect the macro-turbulence 
and have 



a Q = 2cr () a/tt {Q aa * T stsa ] (0) 



(B.l) 



I(x,y) 



(A. 11) 



with T sxsr = To * H and % is the normalized rotation kernel 
( B. 14[ >. If we assume that To is Gaussian with dispersion /?o, the 
associativity of the convolution product leads to 



fl O = 2ff- o -0r[^ <r( *K] (0) 



Griffs,** 

x I duu ®j(x, y, v - xV sin i+) . 

xj — OO 

The inner integral gives 

X+oo 
dvu 2 ® R (u- xVsmi*) = ( 2 cos 2 6 + x 2 (V sin /*) 2 (A. 12) 
CO 

for the radial broadening, and 

X+oo /-*v sin 1+ 

dvu 2 ® T (u- xV sin i+) = £ 2 sin 2 6» + * 2 ( V sin /*) 2 (A.13) a Q = 4cr | Gcr,(v)R(v) du . 
00 Jo 



/V sin 1* 
V sin 



(B.2) 



with erf — ctq + Pq. Finally, since Qa-,( v ) an ^ ^( y ) are even func- 



tions of v, the expression of «o can be slightly simplified 

V sin i+ 



(B.3) 
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The amplitude ao is thus given by one single integral over a 
finite interval, which only has to be computed once. It can be 
done numerically using, for example, the simple trapezoidal rule 
explained in |Press et al.| ( p^92] >. In the case of quadratic limb- 
darkening, the result can also be expressed explicitly as a com- 
bination of modified Bessel functions and error functions (e.g. 
|Hirano et al.| |2010l Eq. F5). 

B.2. With macro-turbulence 



If the macro-turbulence is taken into account, it is not any- 
more possible to express the amplitude qq as a simple integral 
as in Sect. IB. II This is because the rotational-macroturbulence 
broadening kernel cannot be expressed as a convolution product 
(Gray 2005 | l. But since ao brings only a small correction with 
respect to the Gaussian case (ao remains close to 1), we simplify 
the problem and approximate the line profile of the nonrotating 
star as a single Gaussian with dispersion /3' Q given by 



with the initial conditions 
b(-l) =n, 

V^r(3/4) 

b(-l/2)= v v ' = 2.3962804694711844... , 

1 j *-rJ 



(B.12) 



£(0) = 2 , 

V^T(l/4) 

&(l/2) = = 1.7480383695280799 .... 

The normalized rotation kernel 7?(u) entering in the expression 
of the amplitude ao ( |B.3| l, or (B.6 1, which is associated to a nor- 
malized intensity 



n 

is then 

7nb{qc n ) 



(B.4) 



m = 2 



V sin ;* 



V sin z* 



2\ 



(a„ + l)/2 



(B.13) 



(B.14) 



instead of two Gaussians with dispersion fio « and flo r defined 
by 

fil^ = f cos 2 6 and p\ T = jSg + £ 2 sin 2 6 . (B.5) 



With this simplification, we recover the expression (B.3 1 where 
cr f as to be replaced by <xj defined by cr' 2 = ctq + /3[ 2 , i.e., 



sin i+ 

a =4cr V?r I @Av)K{v 
Jo 



B.3. Rotation kernel 



) dv . 



(B.6) 



We now give the expression of the rotation kernel 7?(o) present 
in the expression of the amplitude ao (B.3 i and (B.6 1, without or 
with macro-turbulence, respectively. First, we consider a simpler 
kernel Hdv) associated to an intensity I a (x, y) of the form 

I a (x,y) = (l-x 2 -y 2 yl 2 . (B.7) 
With u = v/(V sin/*), the rotation kernel (Gray 2005} reads as 

l\-u 2 



V sin 



i r vl -"■ 

in u J Vw 



(I - u 2 - y 2 ) a/2 dy . 



(B.8) 



To simplify the integral, we make the change of variable z 
yl Vl - u 2 . We obtain 



KAv) 



b(a) 
V sin i+ 



1 



V sin i+ 



2\ 



(Q-+D/2 



where 

b{a)= f (\-z 2 ) al2 dz. 



(B.9) 



(B.10) 



In practice, for the quadratic and the nonlinear limb-darkening, 
only the cases a — n/2, n € N are used. The integrals b(a) can 
thus be computed using the recurrence relation 



b(a) 



a 



a + 1 



b(a - 2) 



(B.ll) 



Appendix C: RM signal measured by the iodine cell 
technique 



In this section, we compute an analytical expansion of undine ( |27] i 
modeling the RM signal measured by the iodine cell technique. 
The expansion is made possible if the subplanet line profile Tp\ a 
and that of the nonrotating star To are both Gaussian. In that 
case, if we denote Qp (v) = To, @p p {v) = T p i a (v), and Hiv) as the 
rotation kernel, we have, on the one hand, 



dT s tav(v) 

dv 



dv 



*R 



00. 



and, on the other hand, 



dT sta 
dv 



pi a 



(Pp) 



dv 



*R 



(P P ) 



(CI) 



(C2) 



with j6 2 = /3q +/3 2 . Thus, both the numerator and the denominator 
involve integrals of the form {dQo-ldv) * R. Let us consider the 
case where % = % a with 



b(a) 



V sin L 



1 - 



V sin i-i 



(C3) 



and 7] — (a + l)/2. The expansion in series of (dQo-ldv) * is 
obtained by expanding fta in the vicinity of v = 0. We have 



b(a) 



Vsim*^ k\ 



z 



(-17)* 



V sin U 



2 k 



(C4) 



where (-rj)k = 1 if k = 0, and (-r])(-T] + 1) . . . (-rj + k - 1) 
otherwise, is the Pochhammer symbol. It should be noted that 
if 77 G N, (—t]\+\ - 0, then the sum is finite and the expansion 
exact. In the following, we consider only truncated sums up to 
an order K, i.e., k < K. We then have 



d@g- 

dv 



(5) = 



XV sin /* 
V sin 1 



1 



b(a) 



2na- V sin /* 



Z 



kl 



(C5) 



V sin L 



2k v - v I (v-v) 2 
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We apply the change of variable x = (v — v)/cr, and set X\ 
-(V sin z* + u)/cr and xi = (V sin i± - D)/cr. We obtain 



d@<r 
dv 



(v) 



b(a) 



V2^cr Vsin/* k=() 



(-rf)k 
kl 



(C6) 



J.v, \Vsmi*) 



The parenthesis inside the integral is then expanded which leads 
to 



dQg- 
dv 



(v) 



1 



2k-m 



V sin z* 



V sin 



. Zj k\ 2-\m) 

* k=0 V m=0 x ' 
\m r-x 

'*/ J.v, 



(C7) 



y ,+1 e-- t2/2 ^ 



The expression is easier to handle when the two sums are in- 
verted. For that, we introduce truncated hypergeometric func- 
tions of the form 



F K (a, b,c;x)= > — — — — , 
*-t (c) k kl 



(C8) 



k=0 



and Q m (i]; x) defined by 



and 



i-rf) 



1 1 



(C9) 
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Qlp+limx) = 2x 
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With these notations, we have 



d@a 
dv 



2K 



1 b(a) 
V2^rVsin., ^ 



QmW, 



(T 



V sin i 



:)T 



V sin 



dx 



(Cll) 



where 



jT^e-^dx = (si gn?/ r2 m / 2 r(l + y) . (C12) 



and y(a, x) = f* ? a ~'e~' dt is the lower incomplete gamma func- 
tion. 

The formula 



ggb of <g7 



C.lll gives the expansion of the numerator 



Moreov er, th e denominator, which is the in- 
tegral of the square of T^J \C.\\ can be computed by numerical 
integration using the same expansion. We observe numerically 



that the convergence of the expansion of iodine, using Eq. ( C. 1 1 
is quite fast. In practice K — 4 already gives accurate results. 
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